% 
% %Edited by Katherine Hudson (November 21, 2014)
% %Original Author: unknown
% 
% %To run this program, you must also have the following programs:
% %exponentialfit.m, isodd.m, and SRfit.m
% 
% %Rename your data file data.txt. It's best to have a separate file for each
% %data set. It's pretty awful. Blame the original author. I'm just polishing
% %the turd here. 

function [eq_modulus, instant_modulus, perm] = StressRelaxFit(struct1, data1, diameter, step_size, height, number_of_steps, L, U, Lower, Upper)

%First step: User Define and Input Sample Info
% diameter = 4;
stepsize = step_size;
h0 = height;
nsteps = number_of_steps;

%% initialization of the code

% data=[points, elapsedtime, scantime ,Displacement(mm), load(g)]
data = zeros(size(data1,1),5); % initialize an array with 5 columns and same number of rows as the txt file
data(:,2) = data1(:,1); % column 2 = elapsed time
data(:,4) = data1(:,2) - data1(1,2); % column 4 = displacement (mm)
data(:,5) = data1(:,3); % column 5 = load(g)

ipt = findchangepts(data1(:,2),'MaxNumChanges',nsteps); % automatically find the steps

%save original data for Stressrelaxfit portion
dataStressfit = data;

%Parse data into time, stress, strain
time = data(:,2);
stress=data(:,5);%.*(.001*9.81*4)/(pi*(diameter/1000)^2); %(:,5)
step = [];%step(:,4)= -(data(:,4)./stepsize);
ini_disp = data(1,4);
% data(:,4)=data(:,4)-ini_disp;
strain=(data(:,4))./h0;     %(:,4)
ri = zeros(size(data1(:,2),1),1);%ri=abs(data(:,4)./stepsize);   %rough incrimenting (for chunking)  (:,4)

ri(1:ipt(1)-1) = 0;
for i=1:size(ipt,1)-1
ri(ipt(i):ipt(i+1)) = i ;  
end
ri(ipt(end):end) = i+1;

data = [time, stress, strain, ri]; % time = col1; stress = load(g); strain = displacement/height; ri = row_index (which step does it belong to)

% PVT - this shit makes no sense; but too afraid to get rid of it.
% (FB) this is already matlab default !(so why set it???)
format short

% repeat the last row of data for no reason and fix the index  (FB??)
data = [data; data(end,:)];
data(end,end) = nsteps + 1;

% round (:,4) to the nearest integer, throwing out rows not within .25
data = data(find(abs(round(data(:,4))-data(:,4)) < 0.25),:);
data(:,4) = round(data(:,4));

% PVT - these two lines actually make sense.  
% find mean stress when ri=(:,4)=0
mean_stress = mean(data(find(data(:,4)==0),2));
% normalize stresses to mean stress
data(:,2) = data(:,2) - mean_stress;

% This removes 'noise'.  Essentially any row that is greater than the row
% that comes after is removed.  To only capture data points where the gel
% is relaxing.  I'm not positive that it does what should be done.  It
% should remove the data points that do not continuing to relax, or I
% actually think that all the data points should be left in.  
% remove rows
data = data(find(data(1:end-1,2) < data(2:end,2)),:);

%% processing

%Chunk data for each step

i=1;
index = zeros(nsteps,1);
chunk_data = {}; % initialize a cell array to hold steps
figure
for i = 1:nsteps

    I = find(data(:,4)==i);
    index(i) = I(1);
    index(i+1) = I(end)+1;

    chunk_data{i} = data(I,1:3);
    
    
    subplot(ceil(nsteps/2),2,i)
    plot(chunk_data{i}(:,1),chunk_data{i}(:,2),'.')
    xlabel('Time(s)')
    ylabel('Load (g)')
    title(sprintf('%d step',i));

end


%Ask user if the plots are ok

% answer=input('Are you happy with the plots?  If so, enter 0.  Otherwise, enter the step number you would like to change\n');
% 
% while answer % if b is nonzero
%     % DO NOTHING
% end

disp('On to the Fit!');

format compact

for i = 1:nsteps
%     [model,real,values,fit]=SRfit(name);
    [model,real,values,fit]=SRfit(chunk_data{i});

    Stress(i,1:4)=model;         %[A,B,C,Strain]    A,B,C are Rawly's curve fit values-have to convert
    R2(i,1)=fit;
    numberplots=nsteps;

    if 1==isodd(nsteps)

        numberplots=nsteps+1;

    end

    subplot(numberplots/2,2,i)
    tit = sprintf('%d step',i);
    plot(real(:,1),real(:,2),'o',values(:,1),values(:,2),'rx-');
    legend('Data', 'Aexp(B*x)+C');
end



A=-1*Stress(:,1);

Tau=-(1./Stress(:,2));

B=Stress(:,3); %-Stress(:1)

stress2=[A,Tau,B,Stress(:,4),R2];

tau=stress2(:,2);

eqstress=stress2(:,3);

inststress=(stress2(:,1)-stress2(:,3));

strain=stress2(:,4);


% graph the Eq Stress-Strain and Inst Stress-strain curve and ask user to select linear region

figure

subplot(2,1,1)

plotstrain = -1*strain;

ploteqstress = -1*eqstress;

plotinststress = -1*inststress;

plot(plotstrain,ploteqstress,'o');

subplot(2,1,2)

plot(plotstrain,inststress,'x');



% ans=input('if all data looks good enter 0.  If not enter points you want to exclude (in row vector if more than one point \n[#,#,#]\n )');

% while ans ~= 0
%     tau(ans',:)=[];
%     eqstress(ans',:)=[];
%     inststress(ans',:)=[];
%     strain(ans',:)=[];
%     figure(2);
%     subplot(2,1,1)
%     plot(strain,eqstress,'o');
%     subplot(2,1,2)
%     plot(strain,inststress,'x');
%     ans=input('If figure 2 is correct enter 0.\n if not re enter the points from figure 1 that you would like to remove.\n example[#,#,#]');
% end


%plot(strain,eqstress,'o');

disp('Choose the linear region of the curve');

% L = input('Lower bounds >>>  ');
% L = 1;

% U = input('Upper bounds >>>  ');
% U = 5;
if L > U;

    disp('Error: Lower bound is not lower than Upper');

    disp('Exiting');

    pause(8);

    exit

end



figure

plot(dataStressfit(:,2),dataStressfit(:,5))

disp('Choose section for time constant analysis');

% Lower = input('Lower time >>>  ');
% Lower = 1;

% Upper = input('Upper time >>>  ');
% Upper = 5;
if Lower > Upper;

    disp('Error: Lower time is not lower than Upper time');

    disp('Exiting');

    pause(8);

    exit

end

% elimate unwanted rows

B = dataStressfit(dataStressfit(:,2)<=Upper,:);

B = B(B(:,2)>=Lower,:);



B(:,5) = -1*B(:,5);

B(:,2) = B(:,2)-Lower;



D = min(B);



B(:,5) = B(:,5)-(D(1,5)+.01);

exponentialfit(B(:,2),B(:,5));



x = 0:(U-L);

y = ans(1).* exp(-ans(2) * x);

lambda = ans(2);

plot(B(:,2),B(:,5),'o',x,y,'-')


% Fit a linear regression to the range specified by user

Fstrain = strain(L:U,1);

Feqstress = eqstress(L:U,1);

p=polyfit(Fstrain,Feqstress,1);

f=polyval(p,Fstrain);

Finststress = inststress(L:U,1);

p2=polyfit(Fstrain,Finststress,1);

f2=polyval(p2,Fstrain);

clf;


% Get value for Ha, Inst and calculate permeability for each step

for i=1:length(tau)

    Ha(i)= (eqstress(i)- p(2))/(strain(i));

    Inst(i)=(inststress(i)-p2(2))/(strain(i));

    perm(i)=(((diameter/1000)/2)^2)./(pi^2*Ha(i)*tau(i));

end

perm=perm';

plotInst = -1*Inst;

plotFstrain = -1*Fstrain;

plotf = -1*f;


%Plot the curve with the fit and the permeability data

subplot(3,1,1);

plot(plotstrain,ploteqstress,'o',plotFstrain,plotf);

title('Eq Stress vs Strain');

subplot(3,1,2);

plot(plotstrain,Ha, 'kp', plotstrain, Ha, 'g-');

title ('Equilibrium Modulus vs. strain');

subplot(3,1,3);

plot(plotstrain, plotInst, 'kp',plotstrain, plotInst);

title ('Instantaneous Modulus vs. strain');

permfinal =(((diameter/1000)/2)^2)./(pi^2*p(1)*(1/lambda));

instant = -1*p2(1);

fprintf('Equilibrium Modulus (Pa) %12.2f\n', p(1))

fprintf('Instanteous Modulus (Pa) %12.2f\n', instant)

fprintf('Permeability %17e\n', permfinal)


end

